ohi logo
OHI-Northeast | OHI Science | Citation policy


1 Data Source

Downloaded: December 14, 2018 (emailed to us by Jeffrey Vieser at NMFS)

Description: Records of Bmsy and Fmsy estimates from stock assessments conducted in the greater Northeast region

Time range: 2004 - 2018

Format: Tabular


2 Data wrangling

Load raw data

#stock assessment data shared by NMFS
#northeast region
ne_data <- readxl::read_excel(file.path(dir_anx, '_raw_data/NOAA_NMFS/stock_assessments/NE_Assessment_Records.xlsx'))
#midatlantic region
midatl_data <- readxl::read_excel(file.path(dir_anx, '_raw_data/NOAA_NMFS/stock_assessments/MA_Assessment_Records.xlsx'))

Clean data

#select just the columns we are interested in

data <- ne_data %>%
  bind_rows(midatl_data) %>%
  select(year = Year, 
         stock = Stock, 
         f_fmsy = `F/Fmsy`, 
         b_bmsy = `B/Bmsy`) %>%
  gather(key = "indicator", value = "value", -year, -stock) %>%
  filter(value > 0) #Smooth skate has F/Fmsy value of -9999 so removing that

ggplot(data, aes(x = year, y = value, color = indicator)) +
  geom_hline(yintercept = 1, color = "darkgray") +
  geom_line() +
  facet_wrap(~stock, labeller = labeller(stock = label_wrap_gen(width = 25)), scales = "free_y") +
  theme_bw() +
  theme(strip.text = element_text(size=8))

Save data

write.csv(data, file = "data/nmfs_stock_assessment_data.csv")

Span plot

stock_ids <- data$stock %>% unique()

spanplot_df <- data %>%
  group_by(stock) %>%
  mutate(yr_min = min(year), yr_max = max(year))

span_plot <- ggplot(spanplot_df, aes(x = stock)) +
  #ggtheme_basic +
  theme(panel.border     = element_blank(),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_line(colour = 'grey90'),
        panel.background = element_blank()) +
  geom_linerange(aes(ymin = yr_min, ymax = yr_max, 
                     color = indicator),
                 position = position_dodge(width = .5), alpha = .8) +
  scale_color_manual(values = c('darkred', 'darkblue'),
                       labels=c("B/Bmsy", "F/Fmsy")) +
  coord_flip() +
  labs(x = 'Stock',
       y = 'Year',
       color = 'Indicator')

span_plot

What species have multiple stocks? These ones will give us a tough time when matching up catch to the stocks.

dups <- data %>%
  separate(stock, into = c("species", "location"), sep = " - ") %>%
  select(-indicator, -value) %>%
  distinct() %>%
  group_by(year, species) %>%
  mutate(count = n()) %>%
  ungroup() %>%
  filter(count > 1) %>%
  mutate(stock = paste0(species, " - ", location)) %>%
  select(stock, species, location) %>%
  distinct()

dups
## # A tibble: 21 x 3
##    stock                             species        location               
##    <chr>                             <chr>          <chr>                  
##  1 Atlantic cod - Georges Bank       Atlantic cod   Georges Bank           
##  2 Atlantic cod - Gulf of Maine      Atlantic cod   Gulf of Maine          
##  3 Haddock - Georges Bank            Haddock        Georges Bank           
##  4 Haddock - Gulf of Maine           Haddock        Gulf of Maine          
##  5 Windowpane - Gulf of Maine / Geo… Windowpane     Gulf of Maine / George…
##  6 Windowpane - Southern New Englan… Windowpane     Southern New England /…
##  7 Winter flounder - Gulf of Maine   Winter flound… Gulf of Maine          
##  8 Winter flounder - Southern New E… Winter flound… Southern New England /…
##  9 Yellowtail flounder - Cape Cod /… Yellowtail fl… Cape Cod / Gulf of Mai…
## 10 Yellowtail flounder - Georges Ba… Yellowtail fl… Georges Bank           
## # … with 11 more rows

Looks like there are just 8.

3 Gapfill

Use a linear model to fill years between assessments.

3.0.1 Also we should remove species with limited time series (e.g. Atlantic Mackerel has just 2005 & 2018 B/Bmsy and only 2018 F/fmsy). Gapfilling these causes a lot of weirdness.

We gapfill at B/bmsy not stock score. If we have 2010 and 2012 data that shows b/bmsy of 0.7 and 1.3 respectively. If we gapfill at B/bmsy, the 2011 value becomes 1. This would get a stock score of it. But if we instead gapfill between 2010 anda 2012 stock scores (0.7 and something like 0.9) we would get 0.8…

data_gf <- data %>%
  group_by(stock, indicator) %>%
  mutate(timeseries = n()) %>% #get number of years for each stock/indicator combo. We need at least 2 years
  filter(timeseries > 1) %>%
  complete(year = 2005:2018) %>%
  ungroup() %>%
  mutate(value2 = ifelse(year == 2018 & is.na(value), zoo::na.locf(value), value), #if 2018 is NA, use 2017 (or last reported value)
         value3 = ifelse(year == 2005 & is.na(value2), zoo::na.locf(value2, fromLast = TRUE), value2)) %>% #if 2005 is NA, roll back the closest value from following years
  group_by(stock, indicator) %>%
  mutate(gf_value = zoo::na.approx(value3),
         gf = ifelse(is.na(value), 1, 0)) %>% #if a value is gapfilled, assign it a 1, otherwise 0
  select(-value, -value2, -value3, -timeseries) %>%
  rename(value = gf_value)

3.1 Calculate stock scores

Convert these value to stock scores.

#grab just B/Bmsy data
  nmfs_b_bmsy <- data_gf %>%
    filter(indicator == "b_bmsy") %>%
    select(stock, year, value)

#grab just F/Fmsy data  
  nmfs_f_fmsy <- data_gf %>%
    filter(indicator == "f_fmsy") %>%
    select(stock, year, value)
  ### 
  b_bmsy_underexploit_penalty <- 0.25
  b_bmsy_underexploit_thresh  <- 3.00
  f_fmsy_underfishing_penalty <- 0.25
  f_fmsy_overfishing_thresh   <- 2.00
  
  ### Apply rolling mean to F/Fmsy
  ## Why do we do this? because B is a less sensitive metric (relies of biological processes) and F can fluctuate pretty easily because it is really just a mgmt decision.
  nmfs_f_fmsy <- nmfs_f_fmsy %>%
    arrange(stock, year) %>%
    group_by(stock) %>%
    mutate(f_fmsy_rollmean = zoo::rollmean(value, k = 3, align = 'right', fill = NA)) %>%
    ungroup() %>%
    select(-value) %>%
    rename(value = f_fmsy_rollmean)
  
  stock_status_layers <- nmfs_b_bmsy %>%
    full_join(nmfs_f_fmsy) %>%
    spread(indicator, value) 
  
########################################################.
##### run each fishery through the Kobe plot calcs #####
########################################################.
### * ram_b_bmsy, ram_f_fmsy
  
  
### Function for converting B/Bmsy values into a 0 - 1 score
  rescale_bprime_crit <- function(fish_stat_df,
                                  bmax, bmax_val) {
    
    ###using NOAA's limits here
    overfished_th  <- 0.8
    ### 
    underfished_th <- 1.2
    
    bmax_adj <- (bmax - underfished_th) / (1 - bmax_val) + underfished_th
    ### this is used to create a 'virtual' B/Bmsy max where score drops
    ### to zero.  If bmax_val == 0, this is bmax; if bmax_val > 0, bmax_adj
    ### extends beyond bmax, to create a gradient where bmax_val occurs at bmax.
    
    fish_stat_df <- fish_stat_df %>%
      # group_by(stock) %>% ### grouping by stock will set b_max by max per stock, instead of max overall
      mutate(b_max     = max(b_bmsy, na.rm = TRUE)) %>%
      ungroup() %>%
      mutate(bPrime = NA,
             bPrime = ifelse(b_bmsy < overfished_th,  ### overfished stock
                             b_bmsy / overfished_th,
                             bPrime),
             bPrime = ifelse(b_bmsy >= overfished_th & b_bmsy < underfished_th,
                             1,                       ### appropriately fished stock
                             bPrime),
             bPrime = ifelse(b_bmsy >= underfished_th,
                             (bmax_adj - b_bmsy) / (bmax_adj - underfished_th), ### underfished stock
                             bPrime),
             bPrime = ifelse(bPrime < 0, 0, bPrime))
    return(fish_stat_df)
  }
  
  
  ### Function to create vertical gradient based on distance from
  ### ideal F/Fmsy value to actual F/Fmsy value
  f_gradient <- function(f, over_f, under_f, fmax, fmin_val) {
    x <- ifelse(f < over_f & f > under_f, 1, NA)
    x <- ifelse(f <= under_f, (f * (1 - fmin_val) / under_f + fmin_val), x)
    x <- ifelse(f >= over_f,  (fmax - f) / (fmax - over_f), x)
    x <- ifelse(f > fmax, NA, x)
    return(x)
  }
  
  ### Function to convert F/Fmsy values into 0 - 1 score
  rescale_fprime_crit <- function(fish_stat_df,
                                  fmax, fmin_val) {
    
    ### params - taken from BC but changed Bcrit to 0 instead of 0.4:
    Bcrit <- 0.5; overfished_th <- 0.8
    ### underfishing_th is set to the idea "1/3 for the birds":
    underfishing_th <- 0.66; overfishing_th  <- 1.2
    
    bcritslope = 1 / (overfished_th - Bcrit)
    ### connecting from (Bcrit, 0) to (overfished_th, 1)
    
    fish_stat_df <- fish_stat_df %>%
      mutate(fPrime = ifelse(b_bmsy < overfished_th & f_fmsy < fmax,
                             f_gradient(f_fmsy + (overfished_th - b_bmsy) * bcritslope,
                                        over_f = overfishing_th,
                                        under_f = underfishing_th,
                                        fmax = fmax,
                                        fmin_val = fmin_val),
                             NA),
             fPrime = ifelse(b_bmsy >= overfished_th & f_fmsy < fmax,
                             f_gradient(f_fmsy,
                                        over_f = overfishing_th,
                                        under_f = underfishing_th,
                                        fmax = fmax,
                                        fmin_val = fmin_val),
                             fPrime),
             fPrime = ifelse(is.na(fPrime), 0, fPrime), ### fill zeros everywhere unscored
             fPrime = ifelse(is.na(f_fmsy), NA, fPrime) ### but if no f_fmsy, reset to NA
      )
    return(fish_stat_df)
  }
  
  stock_status_df <- stock_status_layers %>%
    rescale_bprime_crit(bmax     = b_bmsy_underexploit_thresh,
                        bmax_val = b_bmsy_underexploit_penalty) %>%
    rescale_fprime_crit(fmax     = f_fmsy_overfishing_thresh,
                        fmin_val = f_fmsy_underfishing_penalty) %>%
    mutate(x_prod = ifelse(!is.na(fPrime), (fPrime * bPrime), bPrime),
           basis  = case_when(
             !is.na(fPrime) & !is.na(bPrime) ~ 'F_Fmsy, B_Bmsy',
             is.na(fPrime)  & !is.na(bPrime) ~ 'B_Bmsy only',
             is.na(bPrime)  & !is.na(fPrime) ~ 'F_Fmsy only'
           )) %>%
    dplyr::select(year, stock,
                  score = x_prod,
                  basis,
                  bPrime, fPrime,
                  b_bmsy, f_fmsy) 

Take the average scores for the 5 species with sub-stocks. Just doing a group_by with species and averaging scores will do this for all species (but only change the values for these sub-species).

Remove stocks with just F_Fmsy - these stocks have no stock scores since we don’t have a way to get a score from just F/Fmsy. Also roll all values forward to 2017.

stock_scores <- stock_status_df %>%
  filter(basis != "F_Fmsy only") %>%
  group_by(stock, year) %>%
  mutate(score = mean(score, na.rm = T)) %>%
  select(year, stock, score, bPrime, fPrime, b_bmsy, f_fmsy, basis)
 
write.csv(stock_scores, file = "data/stock_scores_nmfs.csv")  

3.2 Visualize

ggplot(stock_scores, aes(x = year, y = score)) +
  geom_line() +
  theme_bw() +
  facet_wrap(~stock)

kobe_plot <- function(sp){
  
 ss_df <- stock_scores %>%
    filter(stock == sp) %>%
    arrange(year) %>%
    mutate(last_bbmsy = last(b_bmsy),
           last_ffmsy = last(f_fmsy),
           last_datayear = last(year)) %>%
   ungroup()

generate_kobe_df <- function(f_fmsy_max = 2.5,
                             b_bmsy_max = 3.0,
                             reso       = 0.01,
                             bmax_val   = 0,
                             fmin_val   = 0,
                             weighting_b = 1) {

  kobe_raw <- data.frame(stock  = 1,
                     f_fmsy = rep(seq(0, f_fmsy_max, reso), each  = round(b_bmsy_max/reso) + 1),
                     b_bmsy = rep(seq(0, b_bmsy_max, reso), times = round(f_fmsy_max/reso) + 1))

  kobe <- kobe_raw %>%
    rescale_bprime_crit(bmax = b_bmsy_underexploit_thresh,
                        bmax_val = bmax_val) %>%
    rescale_fprime_crit(fmax = f_fmsy_overfishing_thresh,
                        fmin_val = fmin_val) %>%
    mutate(x_geom  = (fPrime * bPrime),
           x_arith = (fPrime + bPrime) / 2)

  return(kobe)
}

bbmsy_lim <- max(round(max(ss_df$b_bmsy, na.rm = TRUE) + .1, 1), 3)
ffmsy_lim <- max(round(max(ss_df$f_fmsy, na.rm = TRUE) + .1, 1), 2.5)
  
kobe_df <- generate_kobe_df(f_fmsy_max = ffmsy_lim,
                           b_bmsy_max = bbmsy_lim,
                           bmax_val = .25,
                           fmin_val = .25)


kobe_stock_plot <- ggplot(data = kobe_df, aes(x = b_bmsy, y = f_fmsy)) +
    theme_bw() +
    geom_raster(alpha = .8, aes(fill = x_geom)) +
    scale_fill_distiller(palette = 'RdYlGn', direction = 1) +
    labs(title = as.character(sp),
         x = 'B/Bmsy',
         y = 'F/Fmsy',
         fill = "Stock score") +
    geom_path(data = ss_df, 
              show.legend = FALSE,
              aes(x = b_bmsy, y = f_fmsy, group = sp),
              color = 'grey30') +
    geom_point(data = ss_df, 
               show.legend = FALSE,
              aes(x = last_bbmsy, y = last_ffmsy)) +
    geom_text(data = ss_df %>%
                mutate(year = ifelse(year/5 == round(year/5) | year == last_datayear, year, NA)), 
              aes(x = b_bmsy, y = f_fmsy, label = year), 
              hjust = 0, nudge_x = .05, size = 2)

return(kobe_stock_plot)
}

Apply function to all species in the stock_scores data frame.

sp <- stock_scores %>%
  filter(basis == "F_Fmsy, B_Bmsy")

plots <- lapply(unique(sp$stock), kobe_plot)
cowplot::plot_grid(plotlist = plots[1:4])

cowplot::plot_grid(plotlist = plots[5:8])

cowplot::plot_grid(plotlist = plots[9:12])

cowplot::plot_grid(plotlist = plots[13:16])

cowplot::plot_grid(plotlist = plots[17:20])

cowplot::plot_grid(plotlist = plots[21:24])

cowplot::plot_grid(plotlist = plots[25:28])

cowplot::plot_grid(plotlist = plots[29:32])

cowplot::plot_grid(plotlist = plots[33:35])

write.csv(stock_scores, file = "data/nmfs_stock_scores.csv")